THE MESON PRODUCTION IN PROTON-PROTON 
COLLISIONS IN NEXT-TO-LEADING ORDER AND 
INFRARED RENORMALONS 



A. I. Ahmadov 1,2 * and R. M. Burjaliyev 2 
1 The Abdus Salam International Centre for Theoretical Physics 
Strada Costiera 11, 34014, Trieste, Italy and 
2 Department of Theoretical Physics, Baku State University 
Z. Khalilov Street 23, AZ-H48, Baku, Azerbaijan 

Abstract 

In this article, we investigate the next-to-leading order contribution of the higher- twist Feynman 
diagrams to the large-pr inclusive pion production cross section in proton-proton collisions and 
present the general formulae for the higher-twist differential cross sections in the case of the running 
coupling and frozen coupling approaches. We compared the resummed next-to-leading order higher- 
twist cross sections with the ones obtained in the framework of the frozen coupling approach and 
leading-twist cross section. The structure of infrared renormalon singularities of the higher twist 
subprocess cross section and it's resummed expression (the Borel sum) are found. It is shown that 
the resummed result depends on the choice of the meson wave functions used in the calculations. 
We discuss the phenomenological consequences of possible higher-twist contributions to the meson 
production in proton-proton collisions in next-to-leading order at RHIC. 

PACS numbers: 12.38.-t, 13.60.Le, 13.87.Fh, 14.40.Aq, 
Keywords: higher-twist, pion wave function, infrared renormalons 



* ahmadovazar@yahoo.com 
Typeset by REVTeX 



1 



I. INTRODUCTION 



The large-order behavior of a perturbative expansion in gauge theories is inevitably dom- 
inated by the factorial growth of renormalon diagrams [1-4]. In the case of quantum chro- 
modynamics (QCD), the coefficients of perturbative expansions in the QCD coupling a s can 
increase dramatically even at low orders. This fact, together with the apparent freedom in 
the choice of renormalization scheme and renormalization scales, limits the predictive power 
of perturbative calculations, even in applications involving large momentum transfers, where 
a s is effectively small. 

A number of theoretical approaches have been developed to reorganize the perturbative 
expansions in an effort to improve the predictability of the perturbative QCD (pQCD). For 
example, optimized scale and scheme choices have been proposed, such as the method of 
effective charges (ECH) [5], the principle of minimal sensitivity (PMS) [6], and the Brodsky- 
Lepage-Mackenize (BLM) scale-setting prescription [7] and its generalizations [8-20]. In 
Ref . [4] , the resummation of the formally divergent renormalon series and the parametrization 
of related higher-twist power-suppressed contributions has been given. 

In general, a factorially divergent renormalon series arises when one integrates over the 
logarithmically running coupling a s (k 2 ) in a loop diagram. Such contributions do not occur 
in conformally invariant theories which have a constant coupling. Of course, in the physical 
theory, the QCD coupling does run. 

Among the fundamental predictions of QCD are asymptotic scaling laws for large-angle 
exclusive processes [21-28]. QCD counting rules were formalized in Refs. [22,23]. These re- 
actions probe hadronic constituents at large relative momenta, or equivalently, the hadronic 
wave function at short distances. In particular, the meson wave function determines the 
leading higher-twist contribution to meson production at high p?. 

The hadronic wave functions in terms of quark and gluon degrees of freedom play an im- 
portant role in the quantum chromodynamics predictions for hadronic processes. In the per- 
turbative QCD theory, the hadronic distribution amplitudes and structure functions which 
enter exclusive and inclusive processes via the factorization theorems at high momentum 
transfers can be determined by the hadronic wave functions, and therefore they are the 
underlying links between hadronic phenomena in QCD at large (nonperturbative) and small 
distances (perturbative). If the hadronic wave functions were accurately known, then we 
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could calculate the hadronic distribution amplitude and structure functions for exclusive 
and inclusive processes in QCD. 

The study of hadron production at large transverse momenta in hadronic interactions 
is a valuable testing ground of the perturbative regime of QCD, providing information on 
both the parton distribution functions (PDFs) in the proton, and the parton-to-hadron 
fragmentation functions (FFs)[29]. 

The frozen coupling constant approach can be applied for investigation, not only for ex- 
clusive processes, but also for the calculation of higher-twist contributions to some inclusive 
processes, for example as large -p? meson photoproduction [30], two-jet +meson production 
in the electron-positron annihilation [31]. In the frozen coupling constant approach was used 
Refs. [30,31] for calculation of integrals, such as 

J 1 - X 

According to Ref.[7], in pQCD calculations, the argument Q 2 of the running coupling con- 
stant should be taken equal to the square of the momentum transfer of a hard gluon in a 
corresponding Feynman diagram, in both the renormalization and factorization scale. But 
defined in this way, a s (Q 2 ) suffers from infrared singularities. For example in Ref. [32], 
Q 2 equals to (xi — l)u and —xit, where il, t are the subprocess's Mandelstam invariants. 
Therefore, in the soft regions x 1 — >■ 0, x 2 — > 0, the integrals in (1.1) diverge and we need 
some regularization methods of a s (Q 2 ) in these regions for their calculation. In Ref. [33], 
the authors investigated the phenomenology of infrared renormalons in inclusive processes. 
The dispersive approach has been devised to extend properly modified perturbation theory 
calculations towards the low-energy region [34]. Connections between power corrections for 
the three Deep Inelastic Scattering sum rules have also been explored in Ref. [35] . 

Investigation of the infrared renormalon effects in various inclusive and exclusive pro- 
cesses is one of the most important and interesting problems in the perturbative QCD. 
It is known that infrared renormalons are responsible for factorial growth of coefficients 
in perturbative series for the physical quantities. But, these divergent series can be re- 
summed by means of the Borel transformation [1] and the principal value prescription [36], 
and effects of infrared renormalons can be taken into account by a scale-setting procedure 
® S (Q 2 ) ~> o: s (exp(f (Q 2 ))Q 2 ) at the one-loop order results. Technically, all-order resum- 
mation of infrared renormalons corresponds to the calculation of the one-loop Feynman 
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diagrams with the running coupling constant a s (— k 2 ) at the vertices or, alternatively, to 
calculation of the same diagrams with nonzero gluon mass. Studies of infrared renormalon 
problems have also opened new prospects for evaluation of power-suppressed corrections 
to processes characteristics [37]. Power corrections can also be obtained by means of the 
Landau-pole free expression for the QCD coupling constant. The most simple and elabo- 
rated variant of the dispersive approach, the Shirkov and Solovtsov analytic perturbation 
theory, was formulated in Ref.[38]. 

By taking these points into account, it may be argued that the analysis of the next- 
to-leading order higher-twist effects on the dependence of the pion wave function in pion 
production at proton-proton collisions by the running coupling approach, are significant 
from both theoretical and experimental [39] points of view. 

In this work we apply the running coupling approach [40] in order to compute the effects 
of the infrared renormalons on the meson production in proton-proton collisions in next-to- 
leading order. This approach was also employed previously[41,42] to calculate the inclusive 
meson production in proton-proton and photon-photon collisions. The running coupling 
method in next-to-leading order pion and kaon electromagnetic form factor was computed 
in [43]. 

As we know, power-suppressed contributions to exclusive processes in QCD, which are 
commonly referred to as higher-twist corrections. The higher-twist approximation describes 
the multiple scattering of a parton as power corrections to the leading-twist cross section. 

We will show that higher-twist terms contribute substantially to the inclusive meson cross 
section at moderate transverse momenta. In addition, we shall demonstrate that higher- 
twist reactions necessarily dominate in the kinematic limit where the transverse momentum 
approaches the phase-spase boundary. 

A precise measurement of the inclusive charged pion production cross section at y/s = 
62 .4 GeV and y/s = 200 GeV is important for the proton-proton collisions program at the 
Relativistic Heavy Ion Collider (RHIC) at the Brookhaven National Laboratory. Another 
important aspect of this study is the choice of the meson model wave functions. In this 
respect, the contribution of the higher-twist Feynman diagrams to a pion production cross 
section in proton-proton collisions is been computed by using various pion wave functions. 
Also, higher-twist contributions which are calculated by the running coupling constant and 
frozen coupling constant approaches are been estimated and compared to each other. Within 
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this context, this paper is organized as follows: In SecJTTl we provide formulas for the 
calculation of the contribution of the high twist diagrams. In Sec. II III we present formulas 
and an analysis of the next-to-leading order higher-twist effects on the dependence of the 
pion wave function by the running coupling constant approach. In Sec. IIVI we provide 
formulas for the calculation of the contribution of the leading-twist diagrams. In Sec. |V] 
we give the numerical results for the cross section and discuss the dependence of the cross 
section on the pion wave functions. We present our conclusions in Sec. I VI I 

II. CONTRIBUTION OF THE HIGH TWIST DIAGRAMS 

The higher- twist Feynman diagrams, which describe the subprocess q\ + q 2 — > 7i + (tt~) + 7 
for the pion production in the proton-proton collision are shown in Fig.l. In the higher-twist 
diagrams, the pion of a proton quark is directly observed. Their 1/Q 2 power suppression 
is caused by a hard gluon exchange between pion constituents. The amplitude for this 
subprocess can be found by means of the Brodsky-Lepage formula [27]: 

M(s,i) = / dxi I dx 2 5(l - Xx - x 2 )^ 7T (xi,x 2 ,Q 2 )T H (s,i;xi,x 2 ). (2.1) 
Jo Jo 

In Eq.(2.1), Tjj is the sum of the graphs contributing to the hard-scattering part of the 
subprocess. The hard-scattering part for the subprocess under consideration is q\ + q 2 — > 
(<?i Q2) + 7, where a quark and antiquark form a pseudoscalar, color-singlet state (qiq 2 ). 
Here Q(xi,x 2 ,Q 2 ) is the pion wave function, i.e., the probability amplitude for finding the 
valence q±q 2 Fock state in the meson carry fractions X\ and x 2 , X\ + x 2 = 1. Remarkably, 
this factorization is gauge invariant and only requires that the momentum transfers in Th 
be large compared to the intrinsic mass scales of QCD. Since the distribution amplitude and 
the hard-scattering amplitude are defined without reference to the perturbation theory, the 
factorization is valid to leading order in 1/Q, independent of the convergence of perturbative 
expansions. 

The hard-scattering amplitude Tjj can be calculated in perturbation theory and repre- 
sented as a series in the QCD running coupling constant a s (Q 2 ). The function $ is intrinsi- 
cally nonperturbative, but its evolution can be calculated perturbatively. In our calculation, 
we have neglected the pion and the proton masses. Turning to extracting the contributions 
of the higher-twist subprocesses, there are many kinds of leading-twist subprocesses in pp 
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collisions as the background of the higher-twist subprocess q\ + qi — > n + (or n ) + 7, such 
as q + q 7 + g(g ->• 7r + (7r-)), g + g ->• 7 + g(g 7r + (7r - )), g + g -> 7 + ->■ 7r + (7r - )) 
etc. The contributions from these leading-twist subprocesses strongly depend on some phe- 
nomenological factors, for example, quark and gluon distribution functions in the proton 
and fragmentation functions of various constituents, etc. Most of these factors have not 
been well determined, neither theoretically nor experimentally. Thus they cause very large 
uncertainty in the computation of the cross section of process pp — > n + (or it') + 7 + X. In 
general, the magnitude of this uncertainty is much larger than the sum of all the higher-twist 
contributions, so it is very difficult to extract the higher-twist contributions. 

The production of a hadron at large transverse momentum, p?, in a hadronic collisions 
is conventionally analyzed within the framework of perturbative QCD by convoluting the 
leading-twist 2 — > 2 hard subprocess cross sections with evolved structure and fragmentation 
functions. The most important discriminant of the twist of a perturbative QCD subprocess 
in a hard hadronic collision is the scaling of the inclusive invariant cross section [22-25,44], 

° mv m E ^ AB -+ cx) = ^^r 1 (2 - 2) 

at fixed xt = 2px/y/s and center-of-mass angle As we know in the original parton model 
[45] the power fall-off is simply n = 4 since the 2 -» 2 subprocess amplitude for point-like 
partons is scale invariant, and there is no dimensional parameter as in a conformal theory. 
The detected hadron C can be produced directly in the hard subprocess reaction as in an 
exclusive reaction, then such direct higher-twist processes can give a significant contribution 
since there is no suppression from jet fragmentation at large momentum fraction carried by 
the hadron, z, and the trigger hadron is produced without any waste of energy [44]. 

The Mandelstam invariant variables for subprocesses q\ + q 2 — > it + (it~) + 7 are defined 

as 

s = (pi+P2) 2 , i=(p 1 -p 7T ) 2 , u = (p l -p 1 f. (2.3) 

In our calculation, we have also neglected the quark masses. We have aimed to calculate 
the pion production cross section and to fix the differences due to the use of various pion 
model functions. We have used five different wave functions: the asymptotic wave func- 
tion (ASY), the Chernyak-Zhitnitsky wave function [28,46], the CLEO wave function [47], 
the Braun-Filyanov pion wave functions [48] and the Bakulev-Mikhailov-Stefanis pion wave 
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function [49]. It should be noted that the wave functions of pions also are developed in 
Refs. [50-52] by the Dubna group: 



C 3/2 (2x-l) + ^C? z (2x-l) 



,3/2, 



$cleo{x,hI) = $ asy (x) [C s /2 (2x - 1) + 0.27C 3/2 (2:r - 1) - 0.22C 4 3/2 (2z - 1) 
$bf(x,i4) = ® asy (x) [C 3 j 2 (2x - 1) + O.UC 3 2 /2 (2x - 1) + 0.25C 4 3/2 (2x - 1) 
®bms(x,vI) = ®asy(x) [c 3 /2 (2x - 1) + O.lSSCf 2 \2x - 1) - 0.13C 4 3/2 (2x - 1) 

$asy(x) = V3f n x(l ~X), 



(2.4) 



C 3/2 (2:r - 1) = 1, CY\2x - 1) = ^(5(2x - l) 2 - 1), 



,3/2, 



C 4 3/2 (2a; - 1) = ^(21(2x - l) 4 - U{2x - I) 2 + 1), 
8 

where f n = 0.923GCV is the pion decay constant. Here, we have denoted by x = x±, the 
longitudinal fractional momentum carried by the quark within the meson. Then, x<i = 1 — x 
and x\ — X2 = 2x — l. The pion wave function is symmetric under the replacement x\ — xi -H- 
x 2 - x 1 . 

Several important nonperturbative tools have been developed which allow specific predic- 
tions for the hadronic wave functions directly from theory and experiments. The QCD sum- 
rule technique and lattice gauge theory provide constraints on the moments of the hadronic 
distribution amplitude. However, the correct pion wave function is still an open problem 
in QCD. It is known that the pion wave function can be expanded over the eigenfunctions 
of the one- loop Brodsky-Lepage equation, i.e., in terms of the Gegenbauer polynomials 
{C 3/2 (2x-l)}: 
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$AX,Q 2 ) = ^asy(x) 



1+ £ a n (Q 2 )C 3 J 2 (2x-l) 



n=2A.. 



(2.5) 



The evolution of the wave function on the factorization scale Q 2 is governed by the 
functions a n (Q 2 ), 

a s (Q 2 ^ WA) 



a n (Q ) = a„(/i ) 



a«(A*o) 



(2.6) 



72 50 74 364 
Po 81' Po 405' % 6 - 
In Eq.(2.6), {7„} are anomalous dimensions defined by the expression, 
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In = Cf 



1 - 



n+1 1 



> + l)(n + 2) 



(2.7) 



The constants a n (/Xg) = a n are input parameters that form the shape of the wave functions 
and which can be extracted from experimental data or obtained from the nonperturbative 
QCD computations at the normalization point ji^. The QCD coupling constant a s (Q 2 ) at 
the one-loop approximation is given by the expression, 



a s (Q 2 ) 



(2.8) 



/3oMQ 2 /A 2 )- 

Here, A is the fundamental QCD scale parameter, (3 is the QCD beta function one-loop 
coefficient, respectively, 

fo = n- 2 -n } . 

The cross section for the higher-twist subprocess q\q 2 — > 7r + (7r~)7 is given by the expression 



At,. . „, &w 2 a E C F [D(t,u)\ 
^ *.«) = — S — 



1 1 



(2.9) 



where 



D(t, u) = e\t \ dx 
Jo 



1 — x 



i 

+ e 2 u I dx 
'o 



^(Q^^Q 2 
1 — X 



(2.10) 



Here = {x\ — l)u, and Q\ = —x±t, represent the momentum squared carried by the 
hard gluon in Fig.l, ei(e 2 ) is the charge of qi(q 2 ) and Cp = |. The higher-twist contribution 
to the large- pr pion production cross section in the process pp — > -n + {it~) + 7 + X is [53]: 

rfcr Z" 1 Z" 1 , , „ / s „ , , s da 



«/0 



dx 1 dx 2 G qi / hl (x 1 )Gg a /h a (x2) — -j(qq ->■ vr7)<5(s + 1 + u). (2.11) 



7T rft 



d 3 p dydp\ ' 
s = x±x 2 s, 

t = Xit, 

u = x 2 u, 
t = —m T yfse~ y = — pry/se~ y , 
u = —rriT\fse y = ~PtVsg v , 



(2.12) 



8 



%2 



x 2 u x 2 p T Jse v 

Xi = = — 

x 2 s + t x 2 s - pry/se v 

Xit Xip T y/se~ y 



x\s + u x\s — pT\fse y ^ 
where ttit - is the transverse mass of pion, which is given by 

2 2,2 

m T = m + p T 

Let us first consider the frozen coupling approach. In this approach, we take the four- 
momentum square Q\ 2 of the hard gluon to be equal the pion's transverse momentum square 
Qi,2 = Pt- f n this case, the QCD coupling constant a s in the integral (2.10) does not depend 
on the integration variable. After this substitution calculation of integral (2.10) becomes 
easy. Hence, the effective cross section obtained after substitution of the integral (2.10) into 
the expression (2.9) is referred as the frozen coupling effective cross section. We will denote 
the higher-twist cross section obtained using the frozen coupling constant approximation by 

We have extracted the following higher-twist subprocesses contributing to the two covari- 
ant cross sections in Eq.(2.11) 

da 1 - , da 2 - , da 3 _ da 4 _ / rt ,n\ 

—j-yud y /)) -jj-(flii->r7j, — r (ud->n'y), — r (du->n'y). (2.13) 
at dt dt dt 

By charge conjugation invariance, we also have 

da 1 , - , . da 3 , _ . da 2 , - , . da 4 , _ . . 

— —{ud — > 7r^7j = — —{ud — )■ 7r 7), and — —(aw — )■ 7r^7j = — ^(a-u — > it 7). (2.14) 
(it dt dt dt 

III. THE RUNNING COUPLING APPROACH AND HIGHER-TWIST MECHA- 
NISM 

In this section, we shall calculate the integral (2.10) using the running coupling constant 
method and also discuss the problem of normalization of the higher-twist process cross 
section in the context of the same approach. 

As is seen from (2.10), in general, one has to take into account not only the dependence 
of a(Qf 2 ) 011 t ne scale Q\ 2 i but also an evolution of $(x, Q\ 2 ) with Q\ 2 . The meson wave 
function evolves in accordance with a Bethe-Salpeter-type equation. Therefore, it is worth 
noting that, the renormalization scale (argument of a s ) should be equal to Q\ = (37 — l)u, 
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Q\ = —xit, whereas the factorization scale [Q 2 in $m(x, Q 2 )} is taken independent from x, we 
assume Q 2 = p\. Such a approximation does not considerably change the numerical results, 
but the phenomenon considered in this article (effect of infrared renormalons) becomes 
transparent. The main problem in our investigation is the calculation of the integral in 
(2.10) by the running coupling constant approach. This integral in the framework of the 
running coupling approach takes the form 

The a s (Xfi 2 Ro ) has the infrared singularity at x — > 1, if A = 1 — x or x — > 0, if A = x and 
as a result integral (3.1) diverges (the pole associated with the denominator of the integrand 
is fictitious, because <3>m ~ (1 — x), and therefore, the singularity of the integrand at x — 1 
is caused only by a s (Xfij io )). For the regularization of the integral, we express the running 
coupling at scaling variable a s (X/ij, ) with the aid of the renormalization group equation in 
terms of the fixed one a s (Q 2 ). The renormalization group equation for the running coupling 
a = a s /ir has the form 



where 



d ^^^J-l H XQ 2 )f-^HXQ 2 )f 



/So = 11 - \n S Pi = 102 - ^n f , 



The solution of Eq.(3.2), with the initial condition 



a(A)| A= i = a = a s (Q 2 )/n, 



is [35] 



(3.2) 



a(X) 



a 



1 + a^- In A 

4 



aft fkia(X) j 4ft/A + a(A) 



~i -i 



(3.3) 



4/3 V a 4/3 //3i + « 

This transcendental equation can be solved iteratively by keeping the leading a k ln k X and 
next-to-leading a k ln k ~ l X powers. For A = (1 — x) these terms are given by 

^fr ln[l + \nX/t] 



(3.4) 



1 + lnA/f 4tt/3 [1 + lnA/t] 2 
The first term in Eq.(3.4) is the solution of the renormalization group Eq.(3.2) with 
leading power accuracy, whereas the whole expression (3.4) is the solution of Eq.(3.2) with 
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next-to-leading power accuracy. After substituting Eq.(3.4) into Eq.(2.10) we get: 



D(t, u) = e\i I dx 



eita 8 (-u) / dxi T -— 

Jo 0--x 1 + 



1 a s (\n 2 Ro )<S> M {.x ) p 2 T ) 



1-x 



+ e 2 w / dx- 





a 



■(VIo)^m(^Pt) 
1-x 



P 2 t) 



a 2 {-u)(3 1 f\ ^M^^^Ml + lnA/^; 



+ mA/* 1 ) 6lt °47r/3 J ^ (1 - x)(l + In A/ti) 2 



e 2 MQ; s (- 1) / axi 



,«JH)A r 1 $ M (x,p 2 T ) ln(l + lnA/t 2 ) 
o -"' 1 (l-x)(l + liiA/t 2 ) 62U 4vr/3 J % (1 - x)(l + In A/t 2 ) 2 



e\ta : 



Jo 



1 + ^9 A a n(f^Q 



q 3 (pt) 



"I 7n/A) 



^ /2 (2x-l) 



;i-x)(l + lnA/ti) 



47T/3, 



o Jo 



1 + E^4,.. a n(/^o) 



<*s(p|) 



7n/A) 



Cl ,2 (2x - 1) 



ln(l + lnA/ti) 



;i-x)(l + lnA/^) 2 



e 2 wa: s 



(-t) / rfx 

Jo 



1 +E?4,.. a n(/ X o) 



a 3 (Q 2 ) 



7n/A) 



2 (2x - 1) 



;i-a:)(l + lnA/t 2 ) 



e 2 u- 



4tt/3 ( 



JO 



1 ^asy(x) 

dx 



1 + EJL a «(/^o) 



"s(Pt) 



7n//30 



C^ /2 (2x-l) 



ln(l + hiA/t 2 ) 



(l-x)(l + lnA/t 2 ) ; 



(3.5) 



where ti = 47r/a s (— w)/3 and £ 2 = 47r/a s (— 1)/? . 
It should be note that in [41], we used t 1 ~ t 2 . 

The integral (3.5) is common and, of course, still divergent, but now it is recast into a 
form, which is suitable for calculation. Using the running coupling constant approach, this 
integral may be found as a perturbative series in a s (Q 2 ) 



oo 
n=l 



(3.6) 



The expression coefficients S n can be written as power series in the number of light quark 
flavors or, equivalently, as a series in power of (3$: 



S n — C n p o 



n— 1 



The coefficients C n of this series demonstrate factorial growth C n ~ (n — 1)!, which might 
indicate an infrared renormalon nature of divergences in the integral (3.5) and corresponding 



11 



series (3.6). The procedure for dealing with such ill-defined series is well known: one has to 
perform the Borel transform of the series [15] 



B[D\{u) = J2 



oo 



n=0 



then invert B[D](u) to obtain the resummed expression (the Borel sum) D(i,u). After this 
we can find directly the resummed expression for D(Q 2 ). The change of the variable x to 
z = ln(l — x), as ln(l — x) = In A. Then 

D(iu) = e,ia(-u)t, C dx e J ?jt^M f d x ^^ Ml + Z M + 

[ , ) 1 s[ )h J dX (l-x)(t 1 + z) Glt 47Tf3 J dX (1- X )( tl + Z y + 

e 2 ua s — t U 2 / dx- — e 2 u — / dx ; r- 3.7 

' Jo (l-x)(t 2 + z) 4vr/3 J (l-x){t 2 + zf 1 ; 

The basic theoretical problem is how to define the Borel sum B[D](u) of the integral in 
Eq.(3.7) for the quantities we are interested in. In QCD this problem is usually solved using 
perturbative methods and calculating the corresponding multiloop Feynman diagrams with 
a one-gluon line, dressed by the chains of fermion bubbles. 

For the calculation the expression (3.7) we will apply the integral representations of 
l/(t + z) v and \n(t + z) / (t + z) 2 [54,55]. After this operation, formula (3.7) is simplified and 
we can extract the Borel sum of the perturbative series (3.6) and the corresponding Borel 
transform in dependence from the wave functions of the meson, respectively. Also after such 
manipulations the obtained expression can be used for numerical computations. 

It is convenient to use the following integral representation for l/(t + z) u and \n(t + z)/ (t + 
zf: 

— / e^+^u^du, Rev > (3.8) 
i u ) Jo 



1 



(t + zf r 

and 



l ^ z y = ^ (t+z)u (± -C- \nu)udu (3.9) 



where C ~ 0.577216 is the Euler-Mascheroni constant. 

After inserting Eq.(3.8) and (3.9) into (3.7). then, we obtain 

" x f°° <S> M (x,pl)e-^ +z ^ u dudx 



D(t : u) = eita s (—u)ti 



o jo 



(1-x) 



diaK-u)^ 2 f 1 [°° <S> M (x,p?r)e- (tl+ ^ u dudx 



fl 

Jo Jo 



4ttA) Jo Jo (1 - x) 
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1 poo 



e 2 ua s (-i)t 2 



§M{x,P 2 T )e- {t2+z)u dudx 
e 2 ua 2 s (-t)/3 1 tl f 1 f°° ^ M (x,p^)e-^ +z ^ u dudx 




o Jo 



4vr/3 



II 

Jo Jo 



(1-x) 



(3.10) 



In the case of <& asy (x) for D(i,u), we get 

4 v / 37r/ 7r ei£ 



D(t,u) 



Po 



f 

JO 



due 



— t\U 



due~ tlU 
I 1 



1 — u 2 — u 



4y/3n f n eiPii 
Pi 



e 2 w 



o 

oo 



due 



I 

Jo 



-tiu 



1 — u 1 — u 
due~ t2U ' ' 
1 



(1 - C - In it - lnti)w+ 

4v / 37r/ 7r e 2/ #iw 
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In the case of the &cz(x, Q 2 ) wave function, we find 
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Also, in the case of the & B ms(%, Q 2 ) wave function, we get 
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Equations(3.11)-(3.14) is nothing more than the Borel sum of the perturbative series 
(3.6), and the corresponding Borel transform in the case Q asy {x) is 
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The series (3.6) can be recovered by means of the following formula 
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The Borel transform B[D\{u) has poles on the real u axis at u — 1; 2; 3; 4; 5; 6, which confirms 
our conclusion concerning the infrared renormalon nature of divergences in (3.6). To remove 
them from Eqs.(3.11)-(3.14) we applied the principal value prescription. 

Hence, the effective cross section obtained after substitution of the expressions (3.10-3.14) 
into the expression (2.10) is referred as the running coupling effective cross section. We will 
denote the higher-twist cross section obtained using the running coupling constant approach 
by (£rr s - 



IV. CONTRIBUTION OF THE LEADING-TWIST DIAGRAMS 

Regarding the higher-twist corrections to the pion production cross section, a compar- 
ison of our results with leading-twist contributions is crucial. We take two leading-twist 
subprocesses for the pion production: (1) quark- ant iquark annihilation qq — > gj, in which 
the n meson is indirectly emitted from the gluon, g — >■ n + {it~) and (2) quark-gluon fusion, 
qg 57, with subsequent fragmentation of the final quark into a meson, q — > 7r + (7r~). The 
corresponding cross sections are obtained in 

gq) = l™ E * s {Q 2 ) e f 2 (t. + tV 





u 


ft* 






i 




i 
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= — hiF^i^i- (4 - 2) 



da, x 7regQ;_Ba s ((5 

For the leading-twist contribution, we find 

Sm = = / dxidx2dz ( G <n/hi( x i)G q2 /h 2 (x 2 )D n g (z)-^ ( ^(qq ->■ # 7 ) + 



where 



G qi/hl (x 1 )G g/h2 (x 2 )D« q (z) — -^(qg -> g 7 ) ) 5(s + t + u), ( 1.3) 

s = XiX 2 s, t = , -u = , z = . (4.4) 

Z Z X\X 2 S 

Dg(z) = Dg + (z) = Dg ( z) and D*(z) represents gluon and quark fragmentation functions 
into a meson containing gluon and quark of the same flavor. In the leading-twist subprocess, 
the n meson is indirectly emitted from the gluon and quark with the fractional momentum 
z. The 5 function can be expressed in terms of the parton kinematic variables, and the z 
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integration can then be done. The final form for the cross section is 

dcr f^~f^~f 1 dcr 

^M= E ^r = ^2 dx ^ dx 2 \G qi / hl {x l )G q2 / h2 {x 2 )Dl(z) ■ ——^(qq^ g-i) + 

G qi/hl (x 1 )G g/h2 (x 2 )D^(z) ■ ~^S19 -> 97)^ = 

xiG qi / hl (x 1 )sx 2 Gg/ ha (x 2 ) g(y - %{qg -> 57)1 • ( 4 - 5 ) 

t at / 

V. NUMERICAL RESULTS AND DISCUSSION 

In this section, we discuss the numerical results for next-to-leading order higher-twist 
effects with higher-twist contributions calculated in the context of the running coupling and 
frozen coupling approaches on the dependence of the chosen meson wave functions in the 
process pp — > 7i + (or 7r~)7 + X. In the calculations, we use the asymptotic wave function 
^asy, the Chernyak-Zhitnitsky <3>cz, the CLEO pion wave function [41], the Braun-Filyanov 
pion wave functions [42], and the Bakulev-Mikhailov-Stefanis pion wave functionf 43]. For 
the higher-twist subprocess, we take qi + q 2 —> {qiq 2 ) + 7 and we have extracted the following 
four higher-twist subprocesses contributing to pp — > n + (or 7r~)7 cross sections: ud — > 7r + 7, 
du — > vr + 7, ud — > 7r~7, du — > 7r~7 contributing to cross sections. For the dominant leading- 
twist subprocess for the pion production, we take the quark-antiquark annihilation qq — » #7, 
in which the n meson is indirectly emitted from the gluon and quark-gluon fusion, qg — > qj, 
with subsequent fragmentation of the final quark into a meson, q — > -n + {it~). As an example 
for the quark distribution function inside the proton, the MRST2003c package [56] has been 
used. The higher twist subprocesses probe the meson wave functions over a large range 
of Q 2 squared momentum transfer, carried by the gluon. Therefore, in the diagram given 
in Fig.l we take Q\ = (27 — l)u, Q\ = —x\t , which we have obtained directly from the 
higher- twist subprocesses diagrams. The same Q 2 has been used as an argument of a s (Q 2 ) 
in the calculation of each diagram. 

The results of our numerical calculations are plotted in Figs. 2- 15. First of all, it is very 
interesting to compare the resummed higher-twist cross sections with the ones obtained in the 
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framework of the frozen coupling approach. In Figs. 2-4 we show the dependence of next-to- 
leading order higher-twist cross sections (J^+Ynlo calculated in the context of the running 
coupling constant approach and the ratios R = (£f + T )-£ /£f+ T )° , (^Ynlo/(K+Ylo 
as a function of the pion transverse momentum p? for different pion wave functions at 
y — 0. It is seen that the values of cross sections (J^^+Ynlo-i an d ratios for fixed y and 
y/s depend on the choice of the pion wave function. As seen from Fig. 2 the next-to-leading 
order higher-twist differential cross section is monotically decreasing with an increase in 
the transverse momentum of the pion. In Figs. 3-5, we show the dependence of the ratios 
R = (K^nIo/^Ylo, (^)n s lo/(K?Yl% (^YnIo/^ L J) as a function of the pion 
transverse momentum for different pion wave functions. Here (E h +Ynloi + )lo + ) 
are the higher-twist cross sections calculated in the context of the running coupling method, 
in the framework of the frozen coupling approach and is the leading-twist cross section, 
respectively. As seen from Fig.5, in the region 2 GeV/c < p T < 4 GeV/c next-to-leading 
order higher-twist cross section calculated in the context of the running coupling method is 
suppressed by about 1-2 orders of magnitude relative to the leading-twist cross section, but in 
the region 5 GeV/c < px < 30 GeV/c is comparable with the cross section of leading- twist. 
In Figs. 6-8, we have depicted next-to-leading order higher-twist cross sections (J^+Ynloi 
ratios R = (^Ynlo/(^Ylo, (£f + T )^£ /(£^) as a function of the rapidity y of the 
pion at xfs = 62.4 GeV and p T = 4.9 GeV/c. At xfs = 62.4 GeV and p T = 4.9 GeV/c, the 
pion rapidity lies in the region —2.52 < y < 2.52. 

As seen from Figs. 6-8, next-to-leading order higher-twist cross section and ratios have 
a different distinctive. In the region (0.2 < y < 2.52) the next-to-leading order higher- 
twist cross section {J^+Ynloi * s suppressed by about 1-2 order of magnitude relative to 
the leading-twist cross section. In the region (—2.52 < y < —1.92), the ratio for all 
wave functions increase with an increase of the y rapidity of the pion and has a max- 
imum approximately at the point y = —1.92. Besides that, the ratio decreases with 
an increase in the y rapidity of the pion. As is seen from Figs. 7-8, the ratio R very 
sensitive to the choice of the meson wave functions. But, as seen from Fig. 8, the ra- 
tio (S^ r )^£ G /(S^+ ) for all wave functions has a minimum approximately at the point 
y = -1.92. Also, the distinction between R($ asy (x)) with R(& CLEO (x, Q 2 )), R($ C z(x, Q 2 )), 
R(&bf(x,Q 2 )) and R(Q BMS (x,Q 2 )) have been calculated. For example, in the case of 
y/s = 62.4 GeV, y = 0, the distinction between R(® asy (x)) with R($i(x,Q 2 )) (i=CLEO, 
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CZ, BF, BMS) as a function of the pion transverse momentum is shown in Table HI Thus, 
the distinction between R{$ asy {x)) and R($i(x, Q 2 ))(i = CLEO,CZ,BF) is maximum at 
Pt = 20 GeV/c, with R($bms{ x )) at = 2 GeV/c but the distinction between R($> asy (x)) 
with R($i(x, Q 2 )){i = CLEO, CZ, BF) is minimum at p T = 2 GeV/c, with J R($ BM5 (a;)) at 

= 20 GeV/c and increase with an increase in py. Such a behavior of i? may be explained 
by reducing all moments of the pion model wave functions to those of $> asy (x) for high 
Q 2 . Also, we have calculated the distinction between R(<& asy (x)) with R{^cleo{x,Q 2 )), 
R(&cz(x,Q 2 )), R(§bf(x,Q 2 )) and R($bms(%, Q 2 )) as a function of the rapidity y of the 
pion. For example, in the case of yfs = 62AGeV, pr = 4.9GeV/c the distinction between 
R{$ asy {x)) with R($i(x,Q 2 )) (i=CLEO, CZ, BF, BMS) as a function of the rapidity y of 
the pion is presented in Table El 

We have also carried out comparative calculations in the center-of-mass energy y/s = 
200 GeV. The results of our numerical calculations in the center-of-mass energies \/s = 
200 GeV are plotted in Figs.9-15. Analysis of our calculations at the center-of-mass ener- 
gies yfs = 62.4 GeV and \/s = 200 GeV, show that with the increase in beam energy values 
of the cross sections, ratio R = (S^4 r )^| /0 /(S^ r )° , and contributions of next-to-leading 
order higher-twist to the cross section decrease by about 1-2 order. Therefore the exper- 
imental investigation of higher-twist effects include renormalon effects conveniently in low 
energy. On the other hand, the higher-twist corrections and ratios R a very sensitive to the 
choice of the pion wave function. Analysis of our calculations at the center-of-mass energies 
y/s = 62.4 GeV and \/s = 200 GeV, show that values of the next-to-leading order cross 
sections decrease by about (25 — 30) percent of magnitude relative to the leading-order cross 
sections. Also, the distinction between R($ asy (x)) with R($cleo( x , Q 2 ))i R(&cz(%,Q 2 )), 
R{^bf{x,Q 2 )) and R(§bms( x , Q 2 )) have been calculated. For example, in the case of 
y/s = 200 GeV, y = 0, the distinction between R($ asy (x)) with R(<S>i(x,Q 2 )) (i=CLEO, 
CZ, BF, BMS) as a function of the pion transverse momentum px is shown in Table IIHI 
Thus, the distinction between R(§ asy (x)) with Q 2 )), (i= CLEO, BF,) is maxi- 

mum at px = 35 GeV/c, with R(§bms(%)) at Pt = 10 GeV/c, but the distinction be- 
tween R($ asy (x)) with R(® C z(x,Q 2 )), R($cleo(x,Q 2 )), R($>bf(x,Q 2 )) is minimum at 
Pt = 10 GeV/c, with R($bms( x )) at pt = 95 GeV/c and increase with an increase in 
Pt- Also, we have calculated the distinction between R(Q asy (x)) with R(Q C leo{. x i Q 2 )), 
R{^cz{x,Q 2 )), R{^bf{x,Q 2 )) and R(&bms( x , Q 2 ))as a function of the rapidity y of the 
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pion. For example, in the case of y/s = 200GeV, px = !5.5GeV/c the distinction between 
R($ asy (x)) with R(<&i(x, Q 2 )) as a function of the rapidity y of the pion is presented in 
Table IVl The calculations show that the ratio R($i(x, Q 2 ))/R(^ asy (x)), (i=CLEO, CZ, 
BF, BMS) for all values of the transverse momentum px of the pion identically equivalent 
to ratio r($;(x, Q 2 ))/r(<$? asy (x)). 

In our calculations of the next-to-leading order higher-twist cross section of the process 
the dependence of the transverse momentum of meson appears in the range of (10~ 5 -j- 
10~ 25 )mb/GeV 2 . Therefore, higher-twist cross section obtained in our paper should be 
observable at RHIC 

VI. CONCLUDING REMARKS 

In this work we have calculated the inclusive meson production via higher-twist mech- 
anism and obtained the expressions for the subprocess qq — > M7 cross section for mesons 
with symmetric wave functions. For calculation of the next-to-leading order cross section we 
have applied the running coupling constant method and revealed infrared renormalon poles 
in the cross section expression. Infrared renormalon induced divergences have been regu- 
larized by the means of the principal value prescripton and the resummed expression (the 
Borel sum) for the next-to-leading order higher-twist cross section has been found. In the 
resummed higher-twist cross section differs considerably from that found using the frozen 
coupling approach, some regions. Also we have demonstrated that next-to-leading order 
higher-twist contributions to meson production cross section in the proton-proton collisions 
have important phenomenological consequences. Our investigation enables us to conclude 
that the higher-twist pion production cross section in the proton-proton collisions depends 
on the form of the pion model wave functions and may be used for their study. Analysis 
of our calculations shows that the magnitude of next-to-leading order cross sections calcu- 
lated in the running coupling approach in some regions is larger than the leading-twist cross 
sections in 1-2 order. 

Further investigations are needed in order to clarify the role of high twist effects in 
this process. We have demonstrated that the resummed result depends on the pion model 
wave functions used in calculations. The proton-proton collisions provide us with a new 
opportunity to probe a proton's internal structure. In particular, meson production in 
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proton-proton collisions takes into account infrared renormalon effects: this opens a window 
toward new types of parton distributions which can not be measured by the deep inelastic 
lepton-proton scatterings. Finally, we discuss the phenomenological consequences of possible 
higher-twist contributions to meson production in proton-proton collisions in next-to-leading 
order at RHIC. Future RHIC measurements will provide further tests of the dynamics of 
large px hadron production beyond leading twist. 
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PT,GeV/c 


R($cleo(x,Q 2 )) 


R($cz(x,Q 2 )) 


R(<3>bf(x,Q 2 )) 


R(^bms{x,Q 2 )) 




K{<2asy[,X)) 




n{<2asy[,X)) 


2 


0.391 


0.212 


0.338 


7.628 


6 


1.236 


0.387 


1.03 


1.632 


20 


5.169 


4.226 


4.326 


2.268 



TABLE I: The distinction between R(<f> asy (x)) with R(<S>i(x,Q 2 )) (i=CLEO, CZ, BF, BMS) at 
cm. energy \/s = 62.4 GeV. 



y 


R(<5>cleo(x,Q 2 )) 


R(<3>cz(x,Q 2 )) 


R(<3>bf(x,Q 2 )) 


R(<5>bms(x,Q 2 )) 


R($ aS y(x)) 


R($ aS y(x)) 


R(^asy(x)) 


R(®asy(x)) 


-2.52 


7.679 


1.398 


7.142 


2.124 


-1.92 


0.241 


0.226 


0.313 


3.486 


0.78 


0.051 


0.637 


4.367 


1.849 



TABLE II: The distinction between R($ asy (x)) with R(<S>i(x,Q 2 )) (i=CLEO, CZ, BF, BMS) at 
cm. energy \fs = 62.4 GeV and pr = 4.9 GeV/c. 



Pt, GeV/c 


R($>cleo(x,Q 2 )) 

R($>asy{x)) 


R(<3>cz(x,Q 2 )) 

R(<t>asy(x)) 


R(<3> BF (x,Q 2 )) 

R(®asy(x)) 


R(^BM S (x,Q 2 )) 
R(<t>asy(x)) 


10 


0.706 


0.246 


0.438 


1.825 


35 


2.746 


0.824 


2.746 


0.684 


95 


2.021 


0.269 


0,807 


0.342 



TABLE III: The distinction between R(<Z> asy (x)) with R(<S>i(x,Q 2 )) (i=CLEO, CZ, BF, BMS) at 
cm. energy y/s = 200 GeV. 
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y 


r($CLEo(x,Q 2 )) 


r(^cz(x,Q 2 )) 




R(^bms(x,Q 2 )) 


r(<P a sy(X)) 




r(<P a sy(X)) 




-2.52 


3.468 


0.564 


2.121 


4.214 


-1.92 


1.068 


0.204 


0.267 


0.768 


0.78 


0.328 


0.579 


3.322 


3.021 



TABLE IV: The distinction between R(<S> asy (x)) with R(<S>i{x,Q 2 )) (i=CLEO, CZ, BF, BMS) at 
cm. energy yfs = 200 GeV and px = 15.5 GeV/c. 
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FIG. 1: Feynman diagrams for the higher-twist subprocess, q\q2 — > vr + (or tt ) r y. 
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FIG. 2: Next-to-leading order higher-twist 7r + production cross section (X'^+Ynlo as a f unc ti° n of 
the pt transverse momentum of the pion at the cm. energy ^/s = 62.4 GeV. 
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FIG. 3: Ratio R = (£f7)-y£o/( s ^)£o> where L0 and NL0 higher-twist contribution are cal- 
culated for the pion rapidity y = at the cm. energy yfs = 62.4 GeV as a function of the pion 
transverse momentum, pt- 
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FIG. 4: Ratio (£^ r )^£ /(£^ r )£ < Q, where higher-twist contribution are calculated for the pion 
rapidity y = at the cm. energy y/s = 62.4 GeV as a function of the pion transverse momentum, 
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FIG. 5: Ratio (S^4 r )^£ /(S^+ ), as a function of the pr transverse momentum of the pion at the 
cm. energy y/s = 62.4 GeV. 
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FIG. 6: Next-to-leading order higher-twist ir + production cross section (SSOjvio' as a function 
of the y rapidity of the pion at the transverse momentum of the pion pr = 4.9 GeV/c, at the cm. 
energy y/s = 62.4 GeV. 




FIG. 7: Ratio );v£o/(^v+ )lo> as a function of the y rapidity of the pion at the transverse 
momentum of the pion p? = 4.9 GeV/c, at the cm. energy y/s = 62.4 GeV. 
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FIG. 8: Ratio (S^ r )^£ /(S^+ ), as a function of the y rapidity of the pion at the transverse 
momentum of the pion pt = 4.9 GeV/c, at the cm. energy yfs = 62.4 GeV . 
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FIG. 9: Next-to-leading order higher-twist 7r + production cross section (£^f0jvz,o as a f unc ti° n of 
the pt transverse momentum of the pion at the cm. energy yfs = 200 GeV. 
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FIG. 10: Ratio R = (^+)n s lo /{^+)lo^ where L0 and NL0 higher-twist contribution are 
calculated for the pion rapidity y = at the cm. energy y/s = 200 GeV as a function of the pion 
transverse momentum, pt- 
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FIG. 11: Ratio (^^)jv£o/(^^)lO' wn ere higher-twist contribution are calculated for the pion 
rapidity y = at the cm. energy y/s = 200 GeV as a function of the pion transverse momentum, 

PT- 
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FIG. 12: Ratio (S^ r )^£ /(S^+ ), as a function of the px transverse momentum of the pion at the 
c.m.energy yfs = 200 GeV . 




FIG. 13: Next-to-leading order higher-twist tt + production cross section (^^+Ynlo-> as a function 
of the y rapidity of the pion at the transverse momentum of the pion pt = 15.5 GeV/c, at the cm. 
energy a/s = 200 GeV. 
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FIG. 14: Ratio (J^+Y^lo/^J^+Ylo^ as a function of the y rapidity of the pion at the transverse 
momentum of the pion pt = 15.5 GeV/c, at the cm. energy y/s = 200 GeV. 
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FIG. 15: Ratio )^£ /(S^+ ), as a function of the y rapidity of the pion at the transverse 

momentum of the pion pt = 15.5 GeV/c, at the cm. energy y/s = 200 GeV. 



33 



